#include "gauss_legendre.h"

/* Armazena os valores iniciais para o método de Gauss-Legendre */
void start_values_gl(mpf_t a, mpf_t b, mpf_t p, mpf_t t)
{
    /* a0 = 1 */
    mpf_init_set_ui(a, 1);

    /* b0 = 1/sqrt(2) */
    mpf_init_set_ui(b, 2); 
	mpf_sqrt(b, b); 
	mpf_ui_div(b, 1, b);

	/* p0 = 1 */
    mpf_init_set_ui(p, 1);    

	/* t0 = 1/4 */
    mpf_init_set_d(t, 0.25);
}

/* Cálculo de a */
void calc_a_gl(mpf_t a1, mpf_t a0, mpf_t b0)
{
	/* a1 = (a0+b0)/2 */   
	mpf_add(a1, a0, b0); 
	mpf_div_ui(a1, a1, 2); 
}

/* Cálculo de b */
void calc_b_gl(mpf_t b1, mpf_t a0, mpf_t b0)
{
	/* b1 = sqrt(a0*b0) */    
	mpf_mul(b1, a0, b0); 
	mpf_sqrt(b1, b1); 
}

/* Cálculo de t */
void calc_t_gl(mpf_t t1, mpf_t t0, mpf_t a1, mpf_t a0, mpf_t p0)
{
	/* t1 = t0-p0*(a0-a1)^2 */    
	mpf_sub(t1, a0, a1); 
	mpf_mul(t1, t1, t1); 
	mpf_mul(t1, t1, p0); 
	mpf_sub(t1, t0, t1); 
}

/* Cálculo de p */
void calc_p_gl(mpf_t p1, mpf_t p0)
{
	/* p1 = 2*p0 */    
	mpf_mul_ui(p1, p0, 2); 
}


/* Cálculo de pi */
void calc_pi_gl(mpf_t pi, mpf_t a, mpf_t b, mpf_t t)
{
	/* pi = (a1+b1)^2/(4*t1) */
    mpf_add(pi, a, b); 
	mpf_mul(pi, pi, pi); 
	mpf_div(pi, pi, t); 
	mpf_div_ui(pi, pi, 4); 
}

/* Gauss-Legendre sequencial */
void gauss_legendre(mpf_t pi, int *iter)
{
	int i;    
	mpf_t pi_ant;
	mpf_t a[2], b[2], p[2], t[2];

    /* Inicialização do método */
    start_values_gl(a[0], b[0], p[0], t[0]);    
	mpf_inits(a[1], b[1], p[1], t[1], pi_ant, NULL);
	
    for(i=0;;i++)
    {
        /* Calcula 'a1', 'b1', 't1', 'p1' */
        calc_a_gl(a[1], a[0], b[0]);
        calc_b_gl(b[1], a[0], b[0]);
        calc_t_gl(t[1], t[0], a[1], a[0], p[0]);
        calc_p_gl(p[1], p[0]);

		/* Armazena o pi anterior para convergência */
        mpf_set(pi_ant, pi);

        /* Cálculo de pi */
        calc_pi_gl(pi, a[1], b[1], t[1]);
	
		gmp_printf("Pi: %.*Ff\n", 15,pi);

        /* Verifica convergência */
        if(converge(pi_ant, pi)) 
			break;
		
		/* Armazena 'a', 'b', 't' e 'p' para a próxima iteração */
		mpf_set(a[0], a[1]);
		mpf_set(b[0], b[1]);
		mpf_set(t[0], t[1]);
		mpf_set(p[0], p[1]);
    }

	/* Número de iterações */
    *iter = i; 

    mpf_clears(a[0], b[0], p[0], t[0], NULL);
	mpf_clears(a[1], b[1], p[1], t[1], pi_ant, NULL);
}

/* Inicializa a estrutura paralela */
void start_struct_gl(g_legendre *gl)
{
	int i;
	gl->conv = FALSE;

	start_values_gl(gl->a[0], gl->b[0], gl->p[0], gl->t[0]);

    /* Vetor para armazenar todos os valores */
	for(i=1; i<26; i++)
		mpf_inits(gl->a[i], gl->b[i], gl->t[i], gl->p[i], NULL);

	mpf_inits(gl->pi[0], gl->pi[1], NULL);
}

/* Limpa estrutura paralela */
void clear_struct_gl(g_legendre *gl)
{
	int i;

	for(i=0;i<26;i++)
		mpf_clears(gl->a[i], gl->b[i], gl->t[i], gl->p[i], NULL);
	 
	mpf_clears(gl->pi[0], gl->pi[1], NULL);
}

/* Thread para o cálculo de 'a' */
void *thread_calc_a_gl(void *arg)
{
    g_legendre *aux = (g_legendre *) arg;

	int id0 = (aux->i);	
	int id1 = (aux->i)+1;

    calc_a_gl(aux->a[id1], aux->a[id0], aux->b[id0]);

    return NULL;
}

/* Thread para cálculo de b */
void *thread_calc_b_gl(void *arg)
{
    g_legendre *aux = (g_legendre *) arg;

	int id0 = (aux->i);	
	int id1 = (aux->i)+1;

    calc_b_gl(aux->b[id1], aux->a[id0], aux->b[id0]);

    return NULL;
}

/* Thread para cálculo de t */
void *thread_calc_t_gl(void *arg)
{
    g_legendre *aux = (g_legendre *) arg;

	/* Começa a armazenar 't' a partir da segunda iteração */
    if(aux->i < 1) 
		return NULL;

	int id0 = (aux->i)-1;	
	int id1 = (aux->i);

    calc_t_gl(aux->t[id1], aux->t[id0], aux->a[id1], aux->a[id0], aux->p[id0]);

    return NULL;
}

/* Thread para cálculo de p */
void *thread_calc_p_gl(void *arg)
{
    g_legendre *aux = (g_legendre *) arg;

	int id0 = (aux->i);	
	int id1 = (aux->i)+1;

    calc_p_gl(aux->p[id1], aux->p[id0]);

    return NULL;
}

/* Thread para cálculo de Pi */
void *thread_calc_pi_gl(void *arg)
{
    g_legendre *aux = (g_legendre *) arg;
	
	/* Começa a armazenar pi a partir da terceira iteração */ 
    if(aux->i < 2) 
		return NULL;

	int id1 = (aux->i)-2;

	/* Armazena o valor de pi para convergência */
	mpf_set(aux->pi[0], aux->pi[1]);
	
	/* Calcula pi */
    calc_pi_gl(aux->pi[1], aux->a[id1], aux->b[id1], aux->t[id1]);

    /* Se convergiu, a flag de convergência se torna verdadeira */
    if(converge(aux->pi[0], aux->pi[1])) 
		aux->conv = TRUE;

    return NULL;
}

/* Gauss-Legendre paralelo */
void gauss_legendre_paral(mpf_t pi, int *iter)
{
	int i, j;
    pthread_t thread[5];
    g_legendre gl;

    /* Inicializa a struct */
	start_struct_gl(&gl);

    /* Realiza as iterações até atingir a precisão */
	for(i=0; !gl.conv ; i++)
    {
        gl.i = i;

        /* Uma thread é criada para cada variável */
        pthread_create(&thread[0], NULL, thread_calc_a_gl, &gl);
        pthread_create(&thread[1], NULL, thread_calc_b_gl, &gl);
        pthread_create(&thread[2], NULL, thread_calc_t_gl, &gl);
        pthread_create(&thread[3], NULL, thread_calc_p_gl, &gl);
        pthread_create(&thread[4], NULL, thread_calc_pi_gl, &gl);

        /* Sincroniza as threads após elas terminarem */
        for(j = 0; j < 5; j++)
            pthread_join(thread[j], NULL);

		gmp_printf("Pi: %.*Ff\n", 15, gl.pi[1]);
    }

	/* Número de iterações */
    *iter = i; 

    clear_struct_gl(&gl);
}
